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1.  Introduction 

The  Department  of  Defense  has  an 
urgent  need  to  both  understand  and  predict 
the  effects  of  high-altitude  (z  >  10  km) 
turbulence  (HAT)  on  systems  that  operate  at 
or  propagate  through  those  altitudes. 
Examples  of  systems  affected  by  HAT 
include  surveillance  aircraft  such  as  the  U-2 
and  the  unmanned  Global  Hawk,  developing 
weapons  systems  such  as  the  Airborne  Laser 
(ABL),  and  developing  communication 
systems  such  as  the  Transformational 
Communication  Satellites  (TSAT).  The 
surveillance  aircraft  are  primarily  affected 
by  mechanical  turbulence,  which  is  the 
fluctuation  of  wind  velocity  and  is  usually 
caused  by  velocity  gradients,  flow  over 
terrain,  or  convection.  High  frequency 
gravity  waves  can  affect  the  aircraft  the 
same  way  and  are  caused  by  the  same 
forcing  mechanisms  as  actual  mechanical 
turbulence  and  are  considered  as  part  of  the 
HAT  forecast  problem.  The  effects  of 
mechanical  turbulence  include  upsetting  the 
autopilot,  degrading  the  quality  of  the 
surveillance  observations,  and  in  the  worse 
case  jeopardizing  the  aircraft  itself.  Laser- 
based  systems  such  as  the  ABL  and  TSAT 
are  affected  by  optical  turbulence,  which  is 
the  fluctuation  of  the  index  of  refraction  and 
is  caused  by  mechanical  turbulence  in  the 
presence  of  temperature  gradients.  For 
directed-energy  weapon  systems  such  as  the 
ABL,  optical  turbulence  causes  the  laser  to 
wander,  spread,  and  scintillate  which  will 
degrade  its  power,  result  in  increased  dwell 
time  to  destroy  a  target,  and  decrease  the 
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effective  range.  For  laser-based 

communication  systems  such  as  TSAT, 
optical  turbulence  can  result  in  significant 
data  dropouts. 

In  order  to  address  the  issue  of 
predicting  high-altitude  turbulence  and  its 
effects  on  Air  Force  systems,  multiple  DoD 
Challenge  projects  have  been  carried  out 
over  the  past  few  years  to  use  high 
resolution  numerical  simulations  to  better 
understand  turbulence  that  occurs  in  the 
stably  stratified  atmosphere  above  the 
planetary  boundary  layer.  Specifically,  two 
approaches  have  been  used.  First,  since  it  is 
generally  acknowledged  that  synoptic  (Ax  - 
100  km)  and  mesoscale  (1  km  <  Ax  <  100 
km)  weather  patterns  are  responsible  for 
generating  high-altitude  turbulence, 
simulations  were  carried  out  nesting  down 
from  synoptic  scale  global  forecast  models 
to  mesoscale  models  to  microscale  models 
(Ax  ~  0(1 00m)).  This  allows  identification 
and  better  understanding  of  the  large  scale 
trigger  mechanisms  that  cause  high 
frequency  waves  and  turbulence  at  high 
altitudes.  Second,  direct  numerical 
simulations  (DNS)  were  conducted  of 
limited  domains  to  isolate  the  onset  and 
evolution  of  turbulence  for  specific  idealized 
situations,  such  as  Kelvin-Helmholtz  (K-H) 
instabilities.  In  Section  2  previous  work 
addressing  this  problem  will  be  summarized. 
Section  3  will  briefly  discuss  the  two  codes 
that  are  currently  being  used  to  carry  out  the 
approaches  discussed  above.  In  section  4 
the  results  from  the  past  year  will  be 
discussed  in  detail.  In  section  5  the 
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implications  of  the  results  presented  here 
will  be  described. 

2.  Previous  Work 

Initial  DNS  of  K-H  instabilities  by 
Weme  and  Fritts  (2000)  visually  depicted 
the  fact  that  optical  turbulence  and 
mechanical  turbulence  do  not  necessarily 
coincide;  mechanical  turbulence  is  large 
when/where  mixing  is  large,  while  optical 
turbulence  can  be  significantly  reduced  by 
mixing.  However,  optical  turbulence  is 
likely  to  be  found  on  the  top  and  bottom 
edges  of  mechanically  mixed  regions  where 
the  isothermal  layer  caused  by  the  mixing 
interacts  with  the  ambient  stably  stratified 
temperature  gradient,  thus  enhancing  the 
temperature  variances  in  those  areas.  Early 
DNS  results  confirmed  observations  of 
balloon  measurements  (Weme  et  al.  2001) 
that  show  that  turbulence  in  the  free 
atmosphere  occurs  in  very  thin  vertical 
layers  with  length  scales  on  the  order  of  1  - 
10  m.  Weme  and  Fritts  (2000)  also 
computed  2nd-order  structure  functions  for 
temperature  and  the  three  component  winds. 
The  results  for  three  wind  components  were 
not  close  to  each  other,  indicating  that  the 
turbulence  was  anisotropic,  non- 
homogeneous  and  thus  non-Kolmogorov. 
This  was  subsequently  confirmed  by  aircraft 
observations  (Wroblewski  et  al.  2007). 

More  recent  simulations  Weme  et  al. 
(2005)  has  expanded  both  the  spatial  and 
temporal  domain  of  the  K-H  instabilities. 
This  work,  leveraged  with  Capability 
Application  Project  (CAP)  grants,  has 
allowed  for  the  simulation  of  the  K-H  under 
different  sub-critical  Richardson  numbers 
(Ri<0.25).  This  has  greatly  increased  the 
understanding  of  the  K-H  process  in  the 
atmosphere.  Particularly  noteworthy  is  the 
striking  difference  in  the  transition  to 
turbulence  under  differing  stratifications. 
However  despite  the  complexity  of  the 
transition  to  turbulence  and  its  dependence 
on  Ri,  there  are  some  simple  laws  that  help 
us  in  our  parameterization  work.  They  are 

a)  the  final  layer  depth  scales  with  Ri'1 2  and 

b)  the  final  Ri  profile  tends  to  0.55. 


Work  with  the  synoptic-scale  to 
microscale  nesting  started  with  relatively 
simple  idealized  cases  such  as  a  jet  stream 
co-located  with  the  tropopause  (Joseph  et  al. 
2004)  using  the  ASU  microscale  model. 
Results  from  this  work  showed  that  with 
fine  spatial  resolution  in  turbulent  regimes, 
one  could  get  good  estimates  of  dissipation. 
Subsequently  the  approach  was  adapted 
(Mahalov  et  al.  2006)  and  applied  using  the 
ASU  microscale  model  nested  within  the 
Weather  Research  and  Forecast  (WRF) 
model.  This  permitted  the  application  of  the 
microscale  code  to  real  cases.  Results  from 
Mahalov  et  al.  (2006)  showed  the  depiction 
of  mountain  waves  off  Greenland  and  the 
depiction  of  thin  layers  of  turbulence  by  the 
microscale  code  that  were  not  detected  by 
the  WRF  model. 

3.  Methodolgy 

The  technical  approach  employed  for 
the  high-resolution  simulations  of 
atmospheric  wind  shears  uses  the  DNS  code 
named  TRIPLE  which  is  an  incompressible 
pseudo-spectral  Navier-Stokes  solver  that 
uses  Fourier  and  Chebyshev  spatial 
discretizations.  It  is  a  mature  production 
code  that  has  served  as  the  workhorse  for  six 
DoD  High  Performance  Computing 
Modernization  Program  (HPCMP)  and  NSF 
Grand  Challenge  projects  and  two  HPCMP 
Capability  Applications  Projects  (CAP) 
Phase  II  projects.  The  numerical  method 
advances  the  Fourier  and/or  Chebyshev 
coefficients  in  spectral  space.  Potential 
truncation  errors  are  removed  using  the 
“tau”-correction  (Weme  1995),  and  vertical 
stacking  of  Chebyshev  domains  is  possible. 
The  method  is  Nth-order  accurate  in  space, 
where  N  is  the  number  of  spectral  modes. 
Depending  on  the  application  (i.e.,  wave¬ 
breaking,  shear,  convection,  etc.),  different 
spatial  representations  are  used.  Possible 
boundary  conditions  include  stress-free,  no¬ 
slip,  fixed-temperature  or  heat-flux, 
radiative,  incident  gravity  waves,  or  any 
combination  thereof.  For  shear-instability 
and  reduced-equation  studies,  greater 
computation  speeds  are  achieved  by 


employing  simpler  boundary  conditions 
permitted  by  these  problems.  For  these 
cases  stress-free  top  and  bottom  boundaries 
permit  Fourier  representations  in  the 
vertical,  and  complicated  inversion 
operations  can  be  avoided. 

Time-stepping  is  conducted  using  the 
mixed  implicit/explicit  3rd-order  Runge- 
Kutta  (RK3)  scheme  of  Spalart  et  al.  (1991), 
which  has  the  same  storage  requirements  of 
most  2nd-order  schemes,  allowing  us  to 
achieve  larger  numerical  domains  than 
would  otherwise  be  possible.  Nonlinear 
terms  are  treated  explicitly,  while  buoyancy 
and  diffusion  are  handled  implicitly. 
Derivatives  (non-linear  terms)  are  computed 
efficiently  in  spectral  (physical)  space. 

Initial  wind-shear  solutions  (i.e.  Weme 
and  Fritts,  2000)  were  computed  in  domains 
of  size  (X,  A/3,  2X)  in  the  (streamwise, 
spanwise,  vertical)  directions,  where  X  is  the 
wavelength  of  the  most  unstable  asymptotic 
linear  eigen  mode  for  the  K-H  instability.  In 
the  initial  period  of  this  Challenge  project 
and  with  additional  resources  from  the  2004 
CAP  II  allocations,  expanded  runs  were 
completed  in  domains  of  size  (4A,  2X,  2X ), 
for  the  Ri  of  0.05,  0.15,  and  0.20  and 
Reynold’s  number  (Re)  of  2500  (Weme  et 
al.  2005).  In  the  past  year  and  again  with 
the  help  of  additional  resources  from  the 
2006  CAP  II  program,  the  0.15  and  0.20 
were  rerun  with  Re  of  2900  and  4000 
respectively.  These  Re  were  chosen  so  that 
all  solutions  required  nearly  the  same 
number  of  spectral  modes,  which  peaked  at 
(nx,  ny,  nz)  =  (3600,  1800,  1800). 

For  the  meso-  to  mircoscale  portion  of 
this  work  the  Advanced  Research  WRF  core 
(Skamarock  et  al.  2005)  was  used.  The 
model  uses  a  time-split  integration  scheme 
where  low-frequency  modes  that  are 
meteorologically  significant  are  integrated 
using  a  RK3  time  integration  scheme  and 
high-frequency  acoustic  modes  are 
integrated  implicitly  in  the  vertical  over 
smaller  time  steps  to  maintain  stability. 
Spatial  discretization  is  done  using  a  C  grid 
staggering  where  normal  velocity  are 
staggered  one-half  grid  length  from  the 
thermodynamic  variables.  Advection  of 


vector  and  scalar  fields  is  in  the  form  of  flux 
divergence,  and  is  performed  using  the  third 
order  Runge-Kutta  time-integration  scheme 
using  a  fifth  and  third  order  accurate  spatial 
discretization. 

The  WRF  model  was  run  for  a  case 
during  the  2005  Terrain-induced  Rotors 
Experiment  (TREx).  The  standard  WRF 
model  was  configured  with  three  nests.  The 
outermost  nest  received  its  initial  and  lateral 
boundary  conditions  from  the  European 
Centre  for  Medium-range  Weather 
Forecasting’s  (ECMWF)  global  forecast 
model  which  provides  data  at  a  horizontal 
resolution  of  25  km  and  91  vertical  levels  up 
to  0.1  hPa.  The  three  nested  WRF  grids  had 
horizontal  resolutions  of  15,  3,  and  1  km  for 
each  succeeding  inner-nest  and  150  vertical 
levels  up  to  10  hPa.  Within  the  third  nest  a 
fourth,  microscale,  nest  is  employed. 
Instead  of  using  the  previous  ASU 
microscale  code,  a  new  code  has  been 
developed,  based  on  the  full  WRF  equations, 
algorithms  and  subroutines.  This  WRF 
compatible  code  runs  in  a  one-way  nesting 
with  WRF  and  allows  for  both  horizontal 
and  vertical  nesting.  Both  upper  vertical 
and  lateral  boundary  conditions  are  relaxed 
to  the  fields  of  the  finest  third  WRF  nest, 
including  the  vertical  velocity.  Following 
microscale  algorithms  from  previous  ASU 
codes,  an  adaptive  staggered  mesh  is 
implemented  in  the  vertical,  with  a  denser 
grid  around  the  tropopause  and  in  the  lower 
stratosphere.  The  new  WRF-like  code  was 
run  with  a  horizontal  resolution  of  333  m 
and  with  720  vertical  levels. 

4.  Results 

a)  K-H  wind  shear  simulations 

Results  from  the  higher  Re  runs  show 
fuller  spectra  and  more  complete  turbulence 
statistics  than  the  previous  runs.  The 
comparisons  of  the  DNS  with  aircraft 
observations  begun  last  year  (Ruggiero  et  al. 
2006)  has  continued  and  has  gone  beyond 
the  validation  of  the  DNS  results.  Analysis 
of  high  frequency  aircraft  measurements  by 
Wroblewski  et  al.  (2007)  has  detailed 


distinct  signatures  in  the  plots  of 
temperature  and  winds  observed  along  the 
aircraft’s  flight  path.  These  plots  have  been 
compared  to  simulated  transects  through  the 
K-H  simulations.  Of  particular  interest  is  to 
compare  the  different  signatures  within  the 
transects  and  identify  particular  signatures 
or  fingerprints  that  are  indicative  of  a 
particular  initial  Ri  number  and  temporal 
point  in  the  turbulence  evolution.  Figure  1 
is  a  comparison  of  the  DNS  and  aircraft 
data.  Both  plots  exhibit  the  classic  cliff- 
ramp  structure  in  the  potential  temperature 
traces.  The  presence  of  the  relatively 
smoothed  ramps  without  steps  indicates  that 
the  aircraft  measurements  occurred  in  the 
latter  stages  of  the  cliff-ramp  evolution. 
Another  comparison  shown  in  Fig.  2  shows 
the  combination  of  cliff-steps  and  cliff 
ramps  that  suggest  either  an  early  time  for 
Ri=0.20  or  late  time  for  Ri=0.15. 
Ultimately  the  Ri  =0.15  scenario  was  chosen 
as  the  best  match  because  of  the  lack  of 
cliffs  in  the  U  wind  field.  The  real 
intriguing  point  of  these  results  is  how  given 
a  set  of  meteorological  variables  (in  this 
case  high-resolution  temperature  and  wind 
data),  one  can  gain  additional  information 
such  as  initial  Ri  and  age  of  the  turbulence 
layer  by  comparing  directly  with  DNS 
results. 

b)  Microscale  nesting 

The  microscale  case  study  highlighted 
here  is  from  the  TREx  Intensive  Observation 
Period  (IOP)  8  that  took  place  on  3 1  March- 
1  April  in  the  Owens  Valley  region  of 
California.  The  primary  focus  of  the 
experiment  was  to  explore  the  structure  and 
evolution  of  atmospheric  rotors  (intense 
low-level  horizontal  vortices  that  form  on 
axis  parallel  to,  and  downstream  of,  a 
mountain  ridge  crest)  as  well  as  associated 
phenomena  in  complex  terrain.  As  such  it 
provides  an  excellent  case  to  study  the 
generation  of  mountains  waves  and  its 
associated  turbulence  at  high  altitudes.  For 
this  case  the  WRF  model  with  the  custom 
WRF  microscale  nest  was  initialized  on  31 
March  2006  at  00  UTC.  Figure  3  shows  a 


cross  section  Richardson’s  number  and  the 
square  of  the  wind  shear  across  the  Owens 
Valley  from  the  microscale  nest  valid  on  1 
April  2006  at  08  UTC.  Clearly  flow  coming 
over  the  over  mountains  from  the  west  is 
causing  increased  shear  and  turbulence  in 
the  valley  itself.  But  what  is  interesting  to 
note  is  that  above  the  approximate  height  of 
the  Sierra  Nevada  mountains  to  the  west  of 
the  valley,  the  atmosphere  is  relatively 
stable  up  to  10-15  km  at  which  point 
turbulence  layers  (Ri  <  0.25)  are  noted 
beginning  in  horizontal  near  the  peak  of  the 
Sierra  Nevada  ridge  and  continuing  down 
wind.  Evidence  of  the  turbulence  layers 
above  10  km  on  the  microscale  nest  is 
confirmed  by  coincident  radiosonde  data 
(not  shown).  The  temperature  profile 
indicates  wave  activity  occurring  above  the 
tropopause  (located  at  approximately  10.5 
km  MSL)  along  with  enhanced  wind  shear. 
Estimations  of  vertical  velocity  based  on  the 
rise  rate  of  the  balloon  show  vertical 
velocities  exceeding  7  m/s  at  several 
altitudes  above  the  tropopause.  These 
values  for  vertical  velocity  would  have  a 
substantial  impact  on  aircraft  operating  at 
these  altitudes. 

5.  Significance  to  DoD 

The  goal  of  this  project  is  to  improve 
prediction  of  high-altitude  turbulence  and 
quantify  its  effects  on  DoD  systems. 
Currently  the  microscale  nesting  work  is 
achieving  the  latter  by  providing  high 
resolution  characterization  of  the 
atmosphere,  including  turbulence,  in  support 
of  testing  of  the  Airborne  Laser  that  has 
begun  this  spring.  Immediately  after  a  laser 
test,  the  WRF  model  with  microscale 
nesting  is  run.  The  results  from  the 
microscale  nest  are  used  to  explicitly 
calculate  turbulence  parameters  such  as  the 
Ri  and  the  refractive  index  structure 
constant,  a  measure  of  optical  turbulence. 
Using  these  variables,  ABL  testers  can 
characterize  the  atmosphere  they  are  lasing 
through  and  predict  future  test  results  and 
operational  capabilities. 


In  the  DNS  results  described  above  it 
was  demonstrated  that  given  a  certain  set  of 
atmospheric  variables,  in  this  case 
temperature  and  winds  along  a  flight  path, 
and  comparing  them  to  DNS  one  can  extract 
additional  information  about  the 
atmosphere,  such  as  Ri  and  the  age  of  the 
turbulence  evolution.  This  leads  to  an 
opportunity  to  infer  the  subgrid-scale 
turbulence  fields  by  comparisons  with  DNS 
results  once  a  set  of  meteorological  variables 
are  obtained  from  an  operational  mesoscale 
model.  This  is  possible  because  of  progress 
in  recent  years  in  applying  Bayesian 
Hierarchical  Modeling  (BHM)  techniques  to 
geophysical  applications  where  the  character 
of  BHM  is  particularly  well  suited  to 
addressing  measurement  and  model 
uncertainty,  while  quantitatively 
incorporating  significant  scientific 
understanding.  An  additional  advantage  of 
the  BHM  approach  is  that  it  produces  a 
probabilistic  forecast  which  is  more  realistic 
for  solving  the  problem  since  mesoscale 
models  cannot  explicitly  resolve  the 
turbulent  flow  and  are  therefore  unable  to 
make  a  deterministic  forecast.  Currently  the 
DNS  results  along  with  the  aircraft  data 
cited  above  and  balloon  measurements  will 
be  used  to  relate  conditional  probabilities  of 
the  occurrence  of  optical  turbulence  based 
on  input  from  operational  mesoscale  models 
to  produce  the  next  generation  of  optical 
turbulence  prediction  for  the  ABL. 
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Figure  1  Transects  of  potential  temperature  and  the  u,  v,  and  w  wind  components  for  (a)  K-H  DNS 
with  Ri  =  0.2  and  non-dimensional  time  t  =  77  and  (b)  aircraft  flights  over  Wales,  U.  K.  6  June  2000. 
For  the  DNS  plots  the  temperature  is  scaled  by  the  product  of  vertical  potential  temperature  gradient 
and  distance  from  the  center  of  the  billow  and  the  winds  are  scaled  by  the  mean  horizontal  wind 
speed. 


